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Fully 3-dimcnsional computations of flow through a long pipe demand a huge number 
of degrees of freedom, making it very expensive to explore parameter space and diffi- 
cult to isolate the structure of the underlying dynamics. We therefore introduce a '2+e' 
dimensional model of pipe flow which is a minimal 3-dimensionalisation of the axisym- 
metric case: only sinusoidal variation in azimuth plus azimuthal shifts are retained, yet 
the same dynamics familiar from experiments are found. In particular the model retains 
the subcritical dynamics of fully resolved pipe flow, capturing realistic localised 'puff '-like 
structures which can decay abruptly after long times, as well as global 'slug' turbulence. 
Relaminarisation statistics of puffs reproduce the memoryless feature of pipe flow and 
indicate the existence of a Reynolds number about which lifetimes diverge rapidly, pro- 
vided that the pipe is sufficiently long. Exponential divergence of the lifetime is prevalent 
in shorter periodic domains. In a short pipe, exact travelling-wave solutions are found 
nearby to flow trajectories on the boundary between laminar and turbulent flow. In a 
long pipe, the attracting state on the laminar-turbulent boundary is a localised structure 
which resembles a smoothened puff. This 'edge' state remains localised even for Reynolds 
numbers where the turbulent state is global. 



1. Introduction 

Laminar flow through a pipe is possible under controlled laboratory conditions up to 



flow rates well beyond those at which turbulence is typically observed. iPfennigerl (|196lh 
achieved laminar flow at Reynolds numbers, Re :— UD/v, as high as 100 000, where U 
is the mean axial speed, D the diameter and v the kinematic viscosity, indicating that, 
rather than the transition to turbulence being via a linear instability, some other mech- 
anism must be responsible. Given an initial disturbance of sufficiently large amplitude, 
self-sustained turbulence is observed for Re of approximately 2000. This turbulent flow 
exhibits distinct spatial structures at different flow rates. For Re up to around 2250 the 
region of turbulence remains localised, with a length of approximately 20 D, and is re- 
ferred to as a 'puff' ( Wvgnanski fc Champagne] 1973 ). At larger flow rates these puffs 



slowly delocalise by splitting into two or more puffs. At much larger Re, of around 2800, 
the disturbances develop into a rapidly expanding active region of turbulence, referred 
to as a 'slug'. No explanation has been offered that predicts such a progression in flow 
regimes and many issues remain unresolved. 

f Present address : Laboratoire d'Hydrodynamique, Ecole Polytechnique, 91128 Palaiseau, 
France 
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The dynamics of perturbations at transitional Reynolds numbers is believed to be 
strongly influenced by a rapidly incre asing number of branches of exact solutions that 
have been found to appear at these fie dFaisst fc Eckhardtll2003l : IWedin fe KerswelfeOQ/j: 



iKerswelll 120051: IPringle fc Kerswelll l2007t IPringle et all bOOSfK At Re < 1750, puffs ar e 



observed to suddenly and unexpectedly decay in experiments ( Peixinho fc Mulhnll2006l ). 
and it has been suggested t hat the turbulent state wanders between these unstable solu- 
tions before relaminarising (jHof et "aLl l2004tlFaTsst fc Eckhardtll2004l) . The same data also 



suggest that the mean lifetime for a puff become s infinite at Re = 1750 i ndicating that 
the puffs become permanent states at this point ( Peixinho fc Mullin 20061 ). This critical 



Re has been reproduced to within 7% using numerical computations which adopted the 
experimental protocol for initiating the puffs and w orked within a long perio dic pipe of 
50 D, so as to realistically capture the puff structure ( Willis fc Kerswelll 2007a ). However, 
experiments using a different way of initiating the pu ffs and designed to capture longer 
puff transients claim that no such critical Re exists (jHof et all 120061 ). 
presented there in a short w 5 D pipe appear to support this conclusion. 



Simulations also 
The obvious 

question is then whether numerically-simulated turbulence which fills a short pipe has 
the same relaminarisation characteristics as localised puff turbulence captured in longer 
numerical domains. A complete statistical study using fully-resolved 3-dimensional com- 
putations across a spectrum of periodic pipe domains remains prohibitively expensive, 
whereas a survey using a realistic model system could provide a clarifying demonstration 
of difference. 

Evidence has also emerged recently that many of the exact solutions known thus far 



region of phase space dSchneider et al. 200761 Kerswell & Tuttv 
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Willis & Kerswell 20081) (see Itano & Tohl (2001): Wang et all ( 


2007 


klViswanathl d2007f) 



for similar observations in channel and plane Couette flow). ISchneider et all (|2007M ) 
have examined the dynamics of flow restricted to lie in this separatrix in a short 5 D 
pipe, finding a t long times a chaotic att ractor apparently centred on a simple travelling 
wave solution (jPringle fc KersweUll2007l ). In such a short pipe, the turbulence naturally 
fills the pipe when triggered and the laminar-turbulent boundary end-state, or 'edge' 
state, is also a global state. In a longer 25D) pipe, however, localised puffs are the 
naturally triggered state at low Re which raises the issue of what the corresponding 'edge' 
state is and how it varies with Re. For example, is it initially localised and does it lose 
localisation at the same Re as the turbulent puff? Again, a realistic model system can 
suggest probable answers to these questions quickly. 

The use of model systems is well established in plane Couette flow, which exhibits the 
same abrupt subcritical transition behaviour as pipe flow. Several approaches have been 
designed to reduce the number of degrees of freedom of th is problem in order dev elop 
more tractable models. The minimal flow unit introduced bv I Jimenez fc Morn ( 1991 ) has 



been useful in identifyi ng the key compone nts that lead to self-sustaining turbul e nce in 



a very small domain (jHamilton et allll995[ ). The model bv lLagha fc Mannevilld (|2007h 



severely truncates the degrees of freedom in the cross-stream direction, but captures 
spanwise and streamwise spatial structures observed in plane Couette flow. Using this 
model, attempts have been made to measure the lifetime of localised turbulence and to 
determine the characteristic structures seen during the relaminarisation process itself. 
Currently, such calculations would be prohibitively expensive for fully 3-dimensional 
models. Severe truncation to only a few Fourier modes in the tilted cross-wise direction 
has also proven useful in det ermining the origin of oblique bands in plane Couette flow 
( Barklev fc TuckermarJl2007t ) . 

The aim of this article is to establish a model system which preserves the rich dynamics 
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of pipe flow but reduces the number of degrees of freedom of the system so considerably 
that the two current issues mentioned above can be probed. The price paid for this 
reduction is, of course, a close quantitative match with fully 3-dimensional pipe flow. 
But this is more than counterbalanced by the ability to isolate and explore, for example, 
a 'puff-like' structure in a much more accessible system. The issues to be addressed within 
the reduced system are as follows, (a) Do the rclaminarisation statistics for turbulent puffs 
differ in character between short and long pipes? In particular, do long- pipe simulations 
indicate a critical Re for sustained puffs whereas short-pipe simulations not? And (6) 
what does the attracting 'edge' state in the laminar-turbulent boundary look like in a 
long pipe? Is it localised like a puff and, if so, does it delocalise at the same Re as a puff? 

Pr evious attempts to find such a reduced mod el have focused o n axisymmetric pipe 
flow ( Patera k. Orszadll981 ) and helical pipe flow ( Landmanlll990al f6l). but in both cases 
the subcritical dynamics of pipe flow is not retained. We briefly revisit these calcu- 
lations to search afresh for evidence of turbulent transients before introducing a new 
2+e-dimensional model which retains the salient features of fully 3-dimensional pipe flow. 
The presentation starts by discuss ing the formulation used for the calculations performed 
throughout this and earlier work ([Willis fc KersweHll2007ol 120081 ). 



2. Formulation 

Given diameter D and fixed mean axial speed U, it is numerically convenient to scale 
lengths by \D and velocities by 2U in the Navier-Stokes equations, leading to 

(d t + u- v)« = -Vp + 4- (i + P) * + irV 2 u, 

Re Re 

where the non-dimensional variable (3 is the fractional pressure gradient, additional to 
the laminar flow, required to maintain a steady U. A Reynolds number Re p , based on 
the applied pressure gradient, is given b y Re v = Re (1 + (3). Our numerical formulation 
is based on the potential formulation of iMarquesI (|l990l ). which is further re-expressed 
to ease numerical solution. An averaging operator is introduced in the axial direction, z, 
which is periodic over a length L = 2n/a, 



i r 



P*0) = T I (•)<**■ ( 2 -2) 



The velocity, u, is then expressed in terms of a potential ip = ip(r,6,z), the axially 
independent flow h = h(r, 9) and a purely axially-dependent potential <fi = </K r 7 ®i z \ 

u = hz + V A (zip) + V A V A (£<j>), (2.3) 

such that P z 4> = 0. Writing the nonlinear terms as b = (u- V) u, the governing equations 
become 

(d t ~ t^V 2 ) V 2 V^ = -(1 - P z ) z ■ V A V A b, (2.4) 

where V 2 := V 2 — d zz . For a boundary condition u — g(9, z), conditions on the potentials 
are 

h = 0, = 0, -d r i/j=g , -V 2 h ct> = g z , (2.5) 
i dei) + d rz cj) = g r , dr Z V 2 h ip - i(9 e V 2 V|0 = Ref-V/\b-d zz r-V/\g 
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(see appendix [Al for details). Note that r ■ V A b — on the boundary unless an internal 
body force is added. Variables are expanded in Fourier modes, 

A(r,8,z) = ^^A km (r) exp(iakz + im9). (2.6) 

As the variables are real, their coefficients satisfy the property Akm = A*_ k _ , where 
* indicates the complex conjugate, and therefore only coefficients with m ^ are kept. 
Numerical truncation in k and m is discussed in following sections. The operator P z 
picks out k = modes and (1 — P z ) retrieves k ^ modes. In addition to the boundary 
conditions (|2.5|) . regularity at the axis imposes symmetries on the Fourier modes across 
the axis. For the potentials, each mode is even(odd) in r if m is even(odd). 

The system for h is simple to solve as it is second order and has two boundary conditions 
on h, one at the boundary and a symmetry condition at the axis. The system for ip and 
<j> is more difficult to invert as it is coupled through the boundary condition. To enable 
numerical solution we reformulate the system for ip and cj> into a set of five equations, 
each second-order in r, and use an influence matrix technique to by-pass the coupled 
boundary condition (see appendix IB1 for details). 

Both finite-difference and Chebyshev expansions have been used in radius. The lat- 
ter is better at low radial resolution, but the former involves only banded matrices (a 
9-point stencil is used), requiring less memory, and is faster for high radial resolutions. 
Time discretisation is second-order using Crank-Nicolson for the diffusion term and an 
Euler predictor step for the non- linear terms. Information from a Crank-Nicolson correc- 
tor step is used to control the timestep size. Nonlinear terms, b, are evaluated using the 
pseudo-spectral method and are dealiased using the |-rule. The code was tested to re- 
produce eigen values about the laminar s tate, eigenvalues about nonlinear travelling wave 
soluti ons from | Wedin fc Kerswelll (|2004h , by direct comparison with a primitive variable 
code (|Kerswell fc Tuttvll2007t l during the r elaminarisation of a perturbed travelling wave, 



and to calculate the turbulent statistics of Eggel s et ali (J1994J) 



3. The absence of turbulence in previous models 



The origi nal calculations fo r axisymmetric pipe flow ([Patera fc Orszaelll98lD and for 
helical flow (|Landmanlll990ar ) found no evidence for turbulence or even long transients 
at Re < 4000. Here we show that this conclusion extends to Re as large as 10 5 and for 
huge initial disturbances suggesting that there are no exact unstable solutions beyond 
Hagen-Poiseuille flow within these dynamical subspaces. 

The axisymmetric model is straightforward to simulate using the numerical algorithm 
described above by time stepping only modes with m = 0. In helical flow, variations in 
and z are reduced to a dependency in the one variable £ = + az, where a is the pitch 
of the helix and periodicity over L = 2~k / a is preserved. In this scenario the flow may be 
expanded u{r, £) = ^ u q (r) exp(ig£), corresponding to taking only modes k — m — > q 
in our formulation. Letting u' be the deviation from the laminar flow, random initial 
disturbances of the form 



u = 



r) 2 (a 2 k 2 + m 2 ) 2 a km exp(iafcz + im9), 



(3.1) 



were applied to the flow (after projection onto the space of solenoidal functions to enforce 
incompressibility) , where the components of a km were random numbers in C s.t. \ak m \ = 
1. 

Figure Q] shows time evolution of random disturbances at Re — 10 000 for a pipe of 
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Figure 1. Decay of axisymmetric (dashed) and helical (solid) perturbations, initially of energy 
up to 40% of the laminar profile (Re = 10, 000 and pipe 2-n D long). 



length 2ir diameters. Initial disturbances were normalised for E' k , /Eo up to 0.4, where 
E' k , Q is the energy of the axially-dependent modes and Eq is the energy of the laminar 
profile; i.e. disturbances of up to 40% of the en ergy of the laminar flow were considered. 
By way of comparison, the turbulent test case of lEggels et oZ| (|l994l ) has only E' k ^ Q /E « 
0.014 at Re = 5300. The number of finite difference points in radius and truncation of the 
Fourier modes was (160, ±256) in (r, z) and (r, £); the large radial resolution was required 
to stably solve for such high initial energies. Axisymmetric disturbances show little sign 
of nonlinear interactions at this Re and decay almost monotonically. Five other sets of 
runs for other random disturbances showed similar behaviour. Helical flow is slightly 
more promising showing occasional moments of growth against a dominant decay. 

Helical flow exhibits a linear instabi lity at relatively small rotation rates about the axis, 
which occurs at longer L for large r Re dMackrodtll976f). Th is supercritical bi furcation and 
the subsequent nonlinear waves (jToploskv fe AkvlasHl988c lLandmanlll9906r) provided an 
excellent test of the helical code. Measuring the rotation rate by i?n = \flRe, where the 
angular velocity fl is in units U /D, there is a supercritical bifurcation to helical waves at 
Rq = —52.43, Re = 4000 and L = 16tt D. Figure [2] shows evolution of five random initial 
disturbances of E' k ^. /Eo up to 0.3 in the presence of rotation; truncation at (80, ±128) 
in (r, £). For Rq = —50 the flow quickly returns the the parabolic profile. At Rq = —100 
and —300, well beyond the linear instability, the flow rapidly returns to a finite- amplitude 
helical wave flow. This supercritical behaviour p ersists for modest r o tation s so that no 
long term turbulent transients can be generated. iBarnes fc Kerswell (l2000h have shown 
that these helical waves themselves undergo a supercritical Hopf bifurcation so that 
solutions cannot obviously be traced back to non-rotating flow. Figure [2] suggests that 
disconnected branches which could lead to subcritical turbulence in rotating helical pipe 
flow are unlikely to exist. 
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Figure 2. Large helical disturbances scaled from 4% to 20% of Eo rapidly settle to parabolic 
flow for Rn — —50 and to laminar helical waves for Rq = —100, —300 (Re = 4000, L — 167tD). 

Our calculations suggest that rotating helical pipe flow follows the supercritical route 
to turbulence via a sequence of supercritical bifurcations, rather than the abrupt subcrit- 
ical behaviour of 3-dimcnsional non-rotating pipe flow. Having seen strong decay at Re 
approximately five times that for which turbulence is observed in the laboratory, it ap- 
pears that neither dynamics restricted to helical or axisymmetric subspaces are relevant 
for the observed transition. 



4. A 2+e-dimensional model 

We now introduce a third model which has high resolution in the cross-stream (radial) 
and streamwise (axial) directions, but only a few modes in the spanwise (azimuth) direc- 
tion. The model was chosen to preserve a high radial resolution as streak features close to 
the wall appear to be important in the self-sustaining process as do detachments from the 
wall during the relaminarisation stages of low Reynolds number turbulence. High axial 
resolution was retained to allow the possibility of localised turbulent structures. This left 
only the azimuthal direction in which to reduce the number of degrees of freedom: only 
Fourier modes m — 0, ±mo were considered, which corresponds to a sinusoidal variation 
in azimuth, an azimuthal shift of the sinusoid and a mean mode. 

4.1. Spatial characteristics 

Localised structures, surprisingly similar to puffs, were captured by the 2+e-dimensional 
model. Figure [3] compares a puff structure from a simulation fully resolved in azimuth (all 
m up to ±24) with a 'puff' from the 2+e-dimensional model (m = —3, 0, 3). The plots are 
of the correct aspect ratio but only half of the computational domain is shown. Puffs from 
the model appear to be similar to resolved puffs in both length and structure, having a 
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Figure 3. Axial vorticity in the (r, «)-plane, 1:1 aspect ratio, flow up the page; only 25 D of a 
periodic 50 D domain shown. From left to right: fully 3-dimensional simulations at Re = 2000, 
2300, 2700 and 2 +e- dimensional simulations at Re = 2600, 3200, 4000. Presented for both 
models are localised 'puffs', the early stages of delocalisation by the generation of a second puff 
downstream, and global slug-turbulence at larger Re. Far right: Energy in separate components 
of the velocity as a function of axial position (units D 2 U 2 ) for the model puff snapshot at 
Re = 2600. 
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Figure 4. Initial trials of the 2 +e- dimensional model for several m-fold rotational symmetries. 



smooth upstream region close to the wall, an active turbulent region, and a dissipative 
region downstream (although note the different Re). Only modest radial resolutions 
were required to observe such structures: a spectral resolution of 35 Chebyshev modes 
was used for the calculations of this section. Several calculations were performed with 
a lower resolution of 25 radial modes, but puff structures tended to elongate, requiring 
a longer pipe and thus offsetting the reduction in computation times. The energy plot 
in figure [3] shows that the while the axial deviation from the mean flow is extended, 
20 D or greater, the roll components are highly localised, extending only 5-101?. Such 
locali sation of the active regi on of the flow has been observed in the full 3-dimensional 
case (jWillis fc Kerswellll2008l) . 



Axial resolution was chosen to approximately match the spectral drop-off in r (ap- 
proximately 4 orders in the magnitude of the spectral coefficients, or 8 orders in the 
power spectrum) and was ±384 for L — 16tt D « 50 diameters. Axial resolution was 
changed proportionally for other L considered in the rest of this section, hence keeping 
the smallest resolved scale fixed. Puffs were found to translate within 2% of U, slightly 
faster than in 3 dimensions, where they travel approximately 10% slower. Also shown in 
figure [3] is that the transition from localised to global turbulence is gradual, as observed 
experimentally. At larger Re the puff becomes delocalised, splitting into two or more 
localised turbulent regions with relatively laminar regions in between. A t much larger Re 
the pr oportion of vigorous turbulence is seen to increase, as recorded bv lGilbrech fc Hale 
( 19651) . 



4.2. Temporal characteristics 

Another important feature captured by the model is that loca lised puffs may survive fo r 
long times before a sudden decay as observed experimentally ( Peixinho Sz. Mullin 20061 ). 



Typical transients for different toq ar e shown in figure 2] which indicates that higher 
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Figure 5. Probability of transient surviving to time T in the mo = 3 model for a pipe of ~ 100 
diameters (L = 32ir D); 100 observations per Re. See the text for how the 'best fit' lines were 
drawn. 



rotational symmetries tend to decay more quickly. As structur es of 3-fold rota tional 
symmetry are t he most frequently observ ed for transitional Re (iDugg l ebv et al. 

n ; 



Schneider efah l2007oJ : iKerswell fc Tuttvl 120071: fwillis fc Kerswelll 1200 



2007; 



m 



3 was 



chosen for analysis of the lifetimes of disturbances. For this mo, azimuthal length scales 
are also comparable to the radial length scale in the model. Figure [5] shows the probability 
distribution function for puff lifetimes based upon 100 runs at each of several Re in a 
pipe ss 100 diameters long (L — 32ttD). Sets of initial puff conditions were generated 
from snaps hots of a long run at a suffi cien tly large Re, similar to the annealing procedure 
adopted in lPeixinho fc Mullinl (|2006l ) and I Willis fc Kerswelll (|2007ah . No dependence on 
the initial condition was observed, however, other than in the very early times of the 
transient. For th e model, times of order 10 4 D/U could be achieved, significantly longer 
than achieved in Willis fc Kerswell (|2007af ) and in less computation time. The log-plot 
shows an exponential distribution indicating a memoryless process. 

Slopes in figure El are based on a median time, r, calculated by the bootstrapping 
method adopted in Willis fc Kerswell ( 20076 '). The shorter slopes represent 95% confi- 
dence intervals for the slope of the data at Re — 3000, and are accurately calculated 
by the bootstrapping method. The best-fit line was found to be poor estimator of the 
mean, being overly sensitive to rare events (e.g. the outlier in figure El for Re = 2900), 
and the goodness-of-fit does not provide the correct error estimate. In the bootstrapping 
method, samples of size N are generated by resampling from the original N observa- 
tions, with equal probability of selecting each. This is repeated 100 000 times and the 
distribution of the means of the samples provides an accurate confidence interval for the 
mean of the original data set. When all data fall within a maximum observation time, 
the confidence intervals generated by this method converge to those predicted by the 
central limit theorem. A 95% confidence interval for r is approximately t± It / \f~N . The 
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Global Turbulence Localised Turbulence 

L(D) a b p Re c 

2tv -2.5 x 1(T 3 6.2 

4tt -4.3 x 1(T 3 11.6 

8tt -6.5 x 10" 3 18.1 5.9 ±2.0 « 5000 

16vr 4.6 ± 1.6 3575 

32tt 4.8 ± 1.6 3450 

Table 1. Best fit parameters for an exponential fit 1/r = exp(a_Re + b) and an algebraic fit 
1/r — (Re c — Re) 13 for different lengths of periodic pipe. 



method is particularly useful when the data is truncated in time and therefore a mean 
cannot be directly calculated. The bootstrapping method easily accommodates such data 
by further resampling when a truncated point is chosen. Any additional error associated 
with the extra resampling is reflected in a wider confidence interval. 

Estimating the median lifetime of a puff from the data is subject to two sources of error: 
the initial transition period during which the flow evolves from the initial condition to 
become a puff, and the presence of a final relaminarisation phase. The former is eliminated 
by considering each of the first observations as an initi al cut-off time and by exam ination 
of the effect of the cut-off on the estimator r. See Willis fc Kerswelll (|20076h for an 
example plot of r with confidence intervals vs. number of observations cut. Removing 
the first few observations eliminates the effect of the transient on r but slightly widens 
the confidence interval. The relaminarisation time error was minimised by identifying a 
threshold 3D energy below which the turbulent flow always relaminarises and applying 
the same value to all runs to indicate the end of the puff lifetime (E' k ^ /E = 0.001). 

The median time, r(D/U), is dependent on the parameter Re. As mentioned above, 
whether r diverges to infinity or not at a finite Re is a matter of ongoing debate. In 



labora tory experiments using a pipe with D — 20mm, L — 785 D. iPeixinho fe Mu llin 
(l2006h found eviden ce that r ~ \/(Re c — Re) with Re c — 1750. In contrast, experiments 
bv lHof et all (j200fih for D = 4mm, L — 7500 D, and using a different method to initiate 
the puff, found that r ~ exp(ci-Re), for some constant c\. Numerical experiments by 
IWillis fe Kerswelll (|2007al f6h using well-resolved puffs in a 50 D periodic pipe, however, 
show lifetimes to be significantly different from the exponential scaling and the simple 
pow er —1 was clearly seen wi th Re c = 1870, overestimating the experimental value of 
1750 Peixinho fc Mullinl (2006) by only 7%. Computational and experimental limitations 
have confined observations of r to O(10 3 — 10 4 ) D/U. 

Figure [6] shows the maximum likelihood estimator for 1/r and 95% confidence inter- 
vals across a s pectrum of pipe lengths. Lifetimes an order of magnitude larger than the 
calculations of I Willis fc KerswelJ (|2007cJ ) were possible in a pipe twice as long, L w 100 
diameters. In short pipes (figure Eh) , where turbulence fills the domain and is therefore 
global, the lifetime appears to follow an exponential scaling 1/t = exp(ai?e + b) with 
values for a and b are given in Table [TJ The 87r D (ss 25 diameter) pipe represents a cross- 
over situation in that the turbulence is only truly pipe-filling for the highest two data 
points whereas it is localised for the four lower Re shown in figure [6j a) . An exponential 
fit through just the two higher Re points, however, seems to fit a steepening trend as the 
pipe lengthens suggested by the 27T D and 4ir D data sets. 

The remaining localised-turbulence points for L = 87r D are better fit by an algebraic 




Figure 6. Sensitivity of the lifetime, r (D/U), of transients to the pipe length. The upper plot 
(a) shows the data for pipe lengths 2n, 4n and 87r (in D) and exponential fits, 1/r = exp(a_Re+fe), 
through all points which correspond to global (slug) turbulence. The 4 leftmost data points for 
8n correspond to localised (puff) turbulence and a 1/r ~ (Re c — Re) 13 fit is shown for these. 
The lower plot (b) shows data where localised (puff) turbulence is present for pipe lengths 8n, 
16tt and 32n ~ 100 diameters long. Here the best fit lines take the form 1/r ~ (Re c — Re) 13 
where the 1/r vs (Re c — Re) plot on log- log scales (tnset) shows these fits as straight lines (the 
same best fit line for L = 87r D is plotted in both (a) and (&)). Points with error bars in Re were 
calculated using a fixed pressure gradient as opposed to fixed mass flux. This data is plotted 
using the time-averaged Re for a 32-7T D pipe and falls precisely on the fixed mass-flux data (as 
for all other points, 100 observations were used for each). 
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expression t = l/(Re c — Re) 13 (see Table 1 for Re c and (3). The fit includes two extra 
points at lower Re for the 8ir D pipe shown in figure [6^6). Doubling to a 16ir D (« 50 
diameter) pipe produces a rapid drop in 1/r and a more clearly defined curvature in the 
data. This length increase is significant because the turbulent 'puffs' in this reduced model 
are approximately 2013 long (see figure [3]). Hence, while such puffs will be significantly 
affected by enforced periodicity over ss 25 diameters, this artificial constraint should be 
substantially relaxed in an ss 50 diameter pipe and almost absent in a pipe of ~ 100 
diameters. By way of confirmation, a 32tt D pipe produces very similar relaminarisation 
data to the 16n D pipe with both being fitted well by the relation r ~ l/(Re c — Re) 13 
based upon similar values for the fitting parameters (see Table 1). Even though the expo- 
nent (3 is relatively poorly constrained, it is clearly different from the v alue of 1 which is 
observed in fully 3 - dimen sional simulations [Willis fe Kerswelll (l2007af ) and experiments 



Peixinho fc Mullin (|200d ). This quantitative discrepancy is undoubtedly an artifact of the 



reduced model and, in fact, is ty pical of other models for different flows (jBottin fc Chate 



1998HLagha fc Manneville1l2007l) . What this reduced model does clearly exhibit, however 



is a qualitative change in relaminarisation behaviour when the preferred localised turbu- 
lent state (a 'puff') is allowed to develop. If the pipe is long enough to accommodate a 
'puff', the presence of a critical Re is suggested whereas when the pipe is shorter than 
a 'puff' so that the turbulence is always global, no such critical threshold is suggested. 
This behavioural change is entirely consistent with the seemingly contradictory results 
obtained recentl y using fully 3-dimensi onal Navier-Stokes calcul ations in short 5 D pipes 
(|Hof et Qi.ll200fih and long 50 D pipes (jWillis fc Kerswelll E)07^h . 



4.2.1. Constant pressure gradient versus constant mass flux 

The reduced model was also used to investigate another open issue: is the form of 
driving, being either a constant imposed pressure-difference between the ends of the 
pipe or a consta nt imposed mass -flux, important for the relaminarisation behaviour of 
turbulent puffs? iHof et all f)2006h use a pressure-drop set-u p (constant mean p r essure 



gradient) in their experiments over a very long pipe whereas iPeixinho fc Mullinl (|2006T ) 
suck fluid through their pipe to produce a constant mass flux. In an averaged sense, the 
two methods of driving the flow are equivalent but instantaneously and locally the flow 
dynamics are different — the mass flux rate can fluctuate for a constant mean pressure- 
gradient flow, and local variations in the pressure gradient cause fluctuations in the total 
pressure-drop for a constant flux flow. Hence puffs evolving in the two situations are 
subtly different, and it is therefore a leading question as to whether they possess the 
same relaminarisation behaviour. 

To answer this, a series of constant pressure-gradient runs were performed in the 32-7T D 
pipe using randomly selected initial conditions from the same long puff run as for the 
constant mass flux runs. Again 100 runs were calculated for each data point at a given 
pressure gradient and then the median lifetime plotted as a function of Re in figure [6)^6). 
Note that a horizontal error bar is plotted which indicates ± one standard deviation 
in the mean Re value. At the largest mean pressure-gradient used, the mean flow rate 
was Re = 3050 with fluctuations having one standard deviation of ±25. This error in 
Re should inversely scale with the length of pipe for flow driven by a pressure head , 
provided that the disturbance remains localised. For the long pipe of IHof et~al. ( 2006? ) 



such variations should be insignificant (1 part in 10 4 ) whereas for a short (numerical) 
periodic pipe very large variations of O(10%) can be expected. The new data points sit 
precisely on the 327rD constant mass-flux curve indicating that, for at least the long pipe 
(w 100 diameters), the precise form of driving is unimportant for the probability of puff 
relaminarisation. 
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4.3. Characteristics of the laminar-turbulent boundary: short pipe 

Given an initially laminar flow, small perturbations decay back to the laminar flow and 
larger perturbations develop into turbulence for sufficiently large Re. This naturally leads 
to the question of what characterises the dividing s et of flows, for which a small pertur- 
bation may lead to cither laminar or turbulent flow. Iltano fc Tohl (|200ll ) used a shooting 
method to find such a boundary in channel flow and discovered that the flow trajectory on 
the boundary settled upon what they thought was a travelling wave solution, but which 
was later identified as a slowly varying part of a periodic orbit ( Toh fe Itano 2003). This 
orbit is stable within the manifold of flows on this laminar-turbulent boundary or 'edge'. A 
similar situation is found in plane Couette flow (|Viswanathll2007t ISchneider et q/.ll2008[ ). 
where a si ngle simple attractor is f ound. Pipe flow exhibits different characteristics, how- 
ever, with ISchneider et all (|20076h finding a chaotic attractor in which tra jectories pass 



nearb y to exact travelling wave solutions (jMellibovskv fe Meseguerll2007t iDuguet et al 



20081) . When the laminar-turbulent boundary dynamics is restricted within certain sym- 



metry subspaces, however, simple attractors do emerge (jDuguet et all 12008). 

Pipe flow is not so different from channel flow and yet they display different dynamics 
on the boundary. It is not difficult to imagine therefore, that the severe truncation of our 
model could also lead to a loss of chaotic behaviour on the boundary. In this section we 
show that the model preserves the chaotic end state for trajectories on the boundary and 
that exact solutions exist. This motivates extension of the results to long computational 
domains, where undirected 3D calculations would be prohibitively expensive. 

Boundary or edge trajectories for the model are shown in figure [71 for L — -k D and 
over a range of Re. While a difference between the edge and developed turbulence may 
be seen by a rapid increase in energy, a clearer measure appears to be (3 related to the 
pressure gradient required to maintain the fixed flux (see equation [2jJ , or equivalently 
the friction. The pressure gradient is less than 10% greater than the laminar value when 
on the edge and is also smooth in time, whereas after a sudden increase to turbulence it 
is as rapidly varying as the energy. 

At the lowest Re shown in figure [7j the energy of the edge is highly variable in time. 
At the next Re = 6500 a period of slow variation is observed. Duguet et all ( 2008h have 
recently demonstrated that the edge can be used to find exact solutions by identifying 
phases where the flow has a relatively simple temporal behaviour. Such an episode is 
marked by a triangle where the instantaneous flow was used as an initial condition for 
a Newton-Krylov code. This converged to the exact solution pictured in figure [8]^ a), 
confirming that travelling wave solutions do exist for the model. The flow exhibits fast 
streaks towards the walls and slow streaks are shed towards the centre. The state found is 
also very similar to that marked by a circle at Re — 10 000 when the trajectory for a while 
appears to settle towards a steady (translating) state. As Re increases the variability on 
the edge surprisingly decreases, more so than can be explained by the increasing viscous 
time. The viscous time D 2 /v scales as Re in our time units, and longer trajectories 
at larger Re have been shown in figure [Jj to compensate. A possible explanation for the 
decreasing variability is that the eigenvalues of the unstable directions from the travelling 
wave s decrease with increasing Re as in plane Couette flow ( Viswanathll2007 ; Wang et al. 



2007), thus enabling closer and longer visitations. 



The variability of the boundary in terms of energy or any other chosen amplitude 
measure also, of course, indicates the range of such attributes for initial conditions which 
will trigger turbulen ce. The fact that it is easier to trigger turbulence as Re increases 
(e.g. iHof et a/.ll2003h is reflected in the decrease in the mean boundary energy as Re 
increases seen in figure[8]J&). Error bars in this figure represent the mean and one standard 
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Figure 7. Edge trajectories for L = n D. Exact travelling wave solution found nearby to the 
triangle at Re = 6500 (shown in figure [8]) has very similar structure to that found at the circle 
for Re = 10000. 
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Figure 8. (a) Cross-sections of an exact travelling-wave solution L — ir D, Re = 6500, wavespeed 
c = 1.56317 at z = 0, \h. Dark regions are slower than the laminar profile and white regions 
faster, (b) Normalised energy E' k ^ /Eo on the boundary (squares) and of turbulence (circles). 
The crosses indicate the normalised total disturbance energy E' / Eo on the boundary. The slope 
for the total disturbance energy (and also E' k u. /Eo ) is —3, implying that the disturbance 
amplitude, A ~ Re~ 15 . At Re = 5300, E' k ^ /E m 0.011 in the reduced model and w 0.014 in 
the fully 3D test case. 
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Figure 9. (a) Snapshots of axial vorticity at two different times from a laminar-turbulent 
boundary trajectory. L — 16n D, Re = 4000. 20 D of ~ 50 D shown. The disturbance is localised 
but with an extended upstream region. ( 6) An exact solution show n over two axial periods, 
L = 1.6941?, discovered using the code of Pringlc & Kcrswell (2007). (c) The exact solution of 
figure [8] found by the Newton- Krylov code of lDuguet et al\ (|200^ i7L = it D. 



deviation above and below this mean but nevertheless, a clean scalin g emerges for the 
disturbance amplitude of A ~ Re~ 15 . Mellibovskv fc Meseguerl (120061) fou nd this scaling 



when considering streamwise perturbations, and Peixinho fc Mullfnl (|2007t ) in laboratory 
experiments with obliquely oriented jets. 

4.4. Characteristics of the laminar-turbulent boundary: long pipe 
The bisection procedure (|ltano fc Tohll200lt [Schneider et all 12007 StlDuguet et aLll2008h 



for isolating the laminar-turbulent boundary can equally be applied to flow within a 
long pipe although the computational demands becomes increasingly intensive. The 2+e- 
dimcnsional model is ideal for a reconnaissance of likely behaviour, including identifica- 
tion of the form the edge state takes. With this motivation, trajectories on the laminar- 
turbulent boundary were found for a pipe ss 50 diameters long using the pressure gradient 
to distinguish between edge and turbulent flow states. Interestingly, the attracting edge 
state which emerges looks exactly like a turbulent puff except at the trailing edge (up- 
stream region): compare figures [3J and [9j a) . In a turbulent puff this is the most energetic 
part possessing a f l uctua tion energy level comparable to homogenised slug turbulence (see 



Willis fc Kerswelll (|2008D ) at higher Re. However, the edge state is noticeably smoother 



even at Re = 4000 (figure O compared to a puff at Re — 2600 (figure [3J which has 
finer scales. The similarities between the two are the strong wall structures slanting into 
the axis at the upstream region, a region where the axial vorticity reaches the axis (the 
trailing edge region), followed by a gradual relaminarisation/decay downstream. Inter- 
estingly it is the upstream region that waxes and wanes (compare the two axial vorticity 
snapshots of figure rather than the passive-looking downstream wake. 

The edge state remains localised just like the turbulent puff up to Re = 3000, but 
surprisingly also remains localised for much higher Re when the puff has given way to 
(global) slug turbulence. To emphasise this, the localised edge state shown in figure [HI at 
Re — 4000 was generated starting from the global disturbance of figure [31 The fact that 
the edge state remains puff-like throughout the puff-to-slug transitional Re range suggests 
that slug turbulence is destabilised puff-turbulence rather than being a separate state 
occupying a different part of phase space. This is certainly consistent with simulations in 
which a puff state smoothly evolves into a slug by slowly expanding upstream as well as 
downstream. A corollary of this, of course, is that a puff and a slug cannot coexist at a 
given Re: to our knowledge there are no reported experimental observations to contradict 
this claim. 

Figure [5] also shows the axial vorticity for a travelling wave with shift-and- reflect sym- 
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metry, found using the method of IPringle fc Kerswelll ( 2007 ) and that for the travelling 
wave shown in figured] (which has no special symmetry). Both have axial vorticity slanted 
from the wall into the central axis reminiscent of the turbulent puff and edge state. A 
recent search for coherent fast-streak states within turbulent puffs has indicated that the 
flow transiently resembles travelling wav e states upstream and downstream of the trail- 
ing edge region (| Willis fc Kerswellll2008f ). Given that the energetic trailing edge region is 
absent in the edge state, there seems an even higher likelihood of seeing coherent states 
there. 

Although the findings in this 2+e-dimensional model are only suggestive of what may 
occur in the fully 3-dimcnsional setting, they are sufficiently interesting to motivate a 
fully 3-dimcnsional long-pipe computation. This is currently underway. 



5. Discussion 

In this paper, we have described the numerical formulation used to simulate transi- 
tional pipe flow (jWillis fc Kerswelfe007ll l2008t IPuguet et aLll2008h . This is ba sed upon 
the po loidal-toroidal potential decomposition of the velocity field discussed by iMarquesl 
(1990), where the difficulty of coupled boundary conditions has been by-passed by influ- 
ence matrix methods here. Reducing the system to five simple second-order equations, 
the method is accurate, relatively simple to implement and computationally efficient. 
This has then been used here to explore dynamical subspaces in the hope of finding a 
reduced system to aid understanding. 

No evidence for turbulent tr ansients has been found in axisymmetric pipe flow, con- 
firming an earlier investigati o n (Patera fc Orszaglll981 ). or in helical pipe flow consistent 
with the work by Landman ( 1990a j]). Rotating pipe flow displays a classic supercriti- 
cal bif ur cation route to turbulenc e (|Mackrodtlll976t iToploskv fc Akvlaslll988t lLandman 
Il990i ftl: iBarnes fc Kerswell l2000h which has no bearing on the non-rotating situation. 
A brief search in rotating helical pipe flow failed to find evidence for any disconnected 
subcritical branches of solution which may have reached back to the non-rotating limit. A 
2+e-dimensional model, which represents a minimal 3-dimcnsionalisation of the axisym- 
metric limit, does, however, possess a subcritical transition scenario and all the important 
spatio-temporal characteristics of fully resolved pipe flow at Re of the same order of mag- 
nitude. Localised disturbances, structurally similar to turbulent puffs, are found in the 
model at low Re (« 2600), which slowly delocalise at intermediate Re (« 3200), and 
rapidly expand into slugs at high Re 4000). Exact unstable travelling wave solutions 
also exist within the model and appear to underpin the dynamics in phase space. 

Within this 2+e-dimensional model, the relaminarisation statistics of the puffs have 
been examined in pipes of varying lengths. For pipes long enough to allow localised 
turbulence to manifest itself, a critical Re is suggested by the data above which the 
'puff' becomes sustained. On the other hand, if the pipe is short so that turbulence fills 
the whole domain, the data is consistent with only transient behaviour, i.e. the turbulence 
always dies eventually. This qualitative change in behaviour as the computational domain 
is varied is consistent wit h the seemingly contradictory results fou nd recently in fully 3- 
dimensional simulations (|Hof et aLll2006t I Willis fc Kerswel]|l2007al ). The need to resolve 
the spatial inhomogeneity of the puff state properly in numerical ex periments is clear. 

Reconciling the con clusions drawn from the experimental data sets (jPeixinho fc Mullin 
2006t iHof eToLll2006h remains a challenge. The existence of both types of r — Re scal- 



ings are not mutually exclusive, however, and resolution of the issue may be related to 
the known sensiti vity of the flow to the exac t structure and amplitude of a perturba- 
tion. Experiments (|Darbvshire fc Mullinlll995l) and simulations (IFaisst fc Eckharddl200i 
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Moehlis et all 120041 ) have indicated that the laminar-turbulent boundary has a fractal- 
like structure where very small changes to the perturbation can completely change the 
fate of the flow. Very small changes to the perturbation can completely change the fate of 
the flow; finite lifetimes can then exist in the presence of an attractor because an initial 
disturbance, apparently large enough to trigger a puff, may actuall y not be within its 



basin of attraction (e.g. Moehlis et al. 2004 ; Mullin fc Peixinho| [2006V In addition, if the 



laminar-turbulent boundary is closely intertwined with the attractor in phase space, a tra- 
jectory can easily be nudged out of the attractor by noise effects, such as pipe roughness, 
temperature changes, pipe misalignment and vibrations. Noise-induced relaminarisations 
of established puffs for Re above Re c do not therefore contradict the existence of a critical 
Re. Conversely, noise could also artificially maintain puff turbulence in the absence of an 
attractor. There is clearly a need for further experimentation. 

It is worth remarking that longer transients in the 100 D pipe (i.e. Re closer to Re c ) 
could, in principle, have been calculated given the computational savings available in the 
model. It was found, however, that the puffs begin to delocalise for Re w 3200 indicating 
that by Re c = 3450, the puff has become unstable to slug-like turbulence. The same 
issue occurs in the real system: puffs delocalise to become slugs for Re = 2250-2500 
( Willis fc Kerswell 2008). No claim has been made in the literature that expanding slugs 



are other than permanently sustained once generated. 

The 2+e-dimensional model has also presented an opportunity to probe the possible 
dynamics on the laminar-turbulent boundary in long pipe flow. Calculations indicate that 
the attracting state in this set is a localised puff-like structure which is smoother and 
less energetic in the trailing edge region than its turbulent puff counterpart. Also intrigu- 
ingly, this end state remains localised way beyond in Re when the puff has delocalised. 
This tends to suggest that the turbulent puff still exists as a solution but has become 
unstable to a slug state. The variability of the flow on the laminar-turbulent boundary 
also highlights the variability in initial conditions which can trigger turbulence. Just fo- 
cusing on the mean energy gives a d isturbance amplitude scaling A ~ Re~ l b consistent 
with some numerical computations ( Mellibovskv fc Meseguer 2006f) in a short pipe and 



laboratory expe riments with a carefully s pecified jet configuration, designed to excite a 
coherent vortex ( Peixinho fc Mullin 20071 ). Clearly, exploring how far this realisation can 



be usefully developed is a promising area for future research. 

In conclusion, we have introduced a model of pipe flow severely truncated in its az- 
imuthal degrees of freedom but otherwise fully resolved in the others. This notwithstand- 
ing, the remaining system captures all of the rich dynamical behaviour observed in pipe 
flow but obtained at a fraction of the computational cost for the full 3-dimensional situ- 
ation. It therefore presents a very accessible arena in which to test ideas and gain some 
insight quickly before deciding to invest a considerable effort in the full 3-dimensional 
system. 
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Appendix A. Boundary conditions 

The Navier-Stokes equation (|2.ip plus the boundary condition u = g(9, z) are equiv- 
alent to (|2.4|) provided that on the boundary 



V A 



Re 



V 2 )u + b 



0. 



(Al) 



where n is its normal (see Marqueslll99d ) . This condition ensures the term in the square 
brackets is equal to a gradient, such as the pressure. If not imposed, this term may be any 
\z where V 2 x = 0. If this is not a gradient, then an unknown body force is introduced. 
For the axisymmetric case \ is constant, the curl of \z is then zero, and therefore the 
condition is redundant. Otherwise, from the diffusion term in (|A lj) . using the properties 



V u = -V A V A u and V A V A (£/) = V(<9 2 /) - iV 2 /, one finds 

7-> T7 , V7 , , _.^-> . 1 r • V A V A V A (iV 2 » + d zz r ■ V A u 

d zz r-VAg, 



r • V A V' z u = f ■ V A V A (iV^V) 

s 1pl - -d 6 4>2 -T 



Or 



which leads to the last condition of (|2.5p . The other conditions simply express u — g(6, z) 
with the gauge freedom tfi = on the boundary. 

The simplified system for the axially averaged flow, h(r, 8), arises because there exists 
a closed circuit c that is not simply-connected running along the axial direction. This 
has an associated condition 



• dl = 0, 



which together with P z on the second curl of the Navier-Stokes equations leads to the 
governing equation for h. The boundary condition h — assumes P z (g z ) — 0, which 
otherwise would correspond to a translating pipe or moving frame. 



Appendix B. Solution for coupled boundary conditions 

Each Fourier mode for tp,(j) is expanded as the superposition 

tp(r) = 4)(r) +aip H (r), (B 1) 

4>(r)=4>(r)+b<j> H (r), 

where the coefficients a and b are scalars. Subscripts k and m have been dropped. The 
barred and superscripted functions solve two distinct systems. Firstly, 

(d t --^V 2 )h =-(l-P z )z- VA VA6, (B2) 

V 2 0i = 2 , 

(dt-^V 2 )^ =z- VA6, (B3) 

with boundary conditions 

2 = = O, -&= 5 „ 9^i = 0, {tgj =gg tZ^O ' (B4) 

where V 2 = (l/r)d r + d rr — m 2 /r 2 — a 2 k 2 and V 2 = (l/r)d r + d rr — m 2 /r 2 . This time- 
dependent system is written in matrix-vector form, according to the time and radial 
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discretisation, then inverted sequentially for 
homogenised system is 



(f>i — ► (j) and ipi — ► ip. The second 



(dt-^V 2 )^ = 0, 



(ft 



7^ 



2 



o, 



with boundary conditions 



= 1, 



lH 

h 



= 0, 8r1>? = l, 



ip H = 
9 r ^ ff = 



if to = 
if m ^ 



(B5) 



(Bl 



(B7) 



As this system has no time-dependent b, solutions with superscript H may be pre- 
computed. The original boundary conditions on <fi and ip are satisfied upon reconstruc- 
tion from the barred and superscripted variables. Two boundary conditions are satis- 
fied by construction for all cases: as 4> and (p H satisfy trivial boundary conditions, then 
<j> = 4> + b(f> H — on the boundary; similarly —<f>\ = g z is satisfied automatically. 

For axisymmetric modes, m = 0, the system is of lower order as ip always appears as 
d r i/j, including in the boundary condition. The condition involving b is not required for 
this case (see appendix fA")) . The simplest solution is to add the boundary condition ip — 
so that we may invert for ip for all modes. When m = the remaining two boundary 
conditions which couple the potentials are satisfied by selecting scalars a and b according 
to the following evaluated on the boundary: 

a = -(d r ijj + ge)/d r ip H , b = -{d rz (f> - g r )/d rz (j) H if to = 0. 



(B8) 



For non-axisymmetric modes, m ^ 0, the condition —d r ip = ge is satisfied automati- 
cally. The last two conditions are satisfied by solving the system for a and b evaluated 
on the boundary, 



-d„il>? 



3 r 



l -d 6 





a 






b 





i dgtp + d rz 
Rer -V A 6 



4>~ 9r 

d zz 



if m ^ 0. (B 9) 



As this only requires the inversion of a 2 x 2 matrix and the H- functions are pre-computed, 
this is an inexpensive way to ensure all boundary conditions are simultaneously satisfied 
to machine precision. 
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